EUROPHYSICS LETTERS 22 May 1997 

Europhys. Lett., (), pp. () 



Charge oscillations in Debye-Hiickel theory 

Benjamin P. Lee 1 and Michael E. Fisher 2 

1 Polymers Division, National Institute of Standards and Technology, Gaitherburg, Mary- 
land 20899, USA 

2 Institute for Physical Science and Technology, University of Maryland, College Park, 
Maryland 20742, USA 



(received ; accepted ) 

PACS. 61.20.Qg- Structure of associated liquids: electrolytes, molten salts, etc. 
PACS. 05.20.— y- Statistical mechanics. 
PACS. 05.70.Jk - Critical point phenomena. 



Abstract. - The recent generalized Debye-Hiickel (GDH) theory is applied to the calculation 
of the charge-charge correlation function, Gzz(r). The resulting expression satisfies both (i) the 
charge neutrality condition and (ii) the Stillinger-Lovett second-moment condition for all T and 
Pn, the overall ion density, and (iii) exhibits charge oscillations for densities above a "Kirkwood 
line" in the (pn,T) plane. This corrects the normally assumed DH charge correlations, and, 
when combined with the GDH analysis of the density correlations, leaves the GDH theory as 
the only complete description of ionic correlation functions, as judged by (i)-(iii), (iv) exact 
low-density (pN,T) variation, and (v) reasonable behavior near criticality. 



A complete theory of ionic fluids, and in particular one which describes the ionic critical 
region, remains an outstanding challenge for statistical physics Q. Recent progress has 
been made at the mean-field level by studies based on the Debye-Hiickel (DH) theory of the 
restricted primitive model (RPM) of electrolytes Eq], supplemented with ion pairing, free-ion 
depletion, and the concomitant dipolar- ionic interactions A satisfactory theory must also 
include a description of the ionic correlation functions. However, the conventionally assumed 
DH ion-ion correlation functions, of which there appear to be three varieties, have several 
shortcomings: the most significant is the absence of a diverging density-density correlation 
length at the DH critical point, which renders the theory totally uninformative as regards 
the relevant order-parameter fluctuations. Further, although the predicted charge-charge 
correlation function, Gzz(?), with the familiar screening decay as e~ K ° r /r, is exact in the 
low-density limit, it violates the Stillinger-Lovett second- moment condition Q. Furthermore, 
there is much evidence indicating that the charge correlations become oscillatory at moderately 
high densities j^, |(| , a phenomenon also missed by the simple screening form. 

Similar difficulties arise with the oft-used mean-spherical approximation (MSA), which 
exhibits, in particular, a complete cancellation of Coulombic effects from the density-density 

Typeset using EURO-L4T E X 



2 



EUROPHYSICS LETTERS 



correlation function, GjvAr(r). In this case some improvement has been found by adding an ad 
hoc term to the assumed direct correlation functions and adjusting it to gain consistency with 
various sum rules; this generalized MSA, or GMSA, yields non-trivial density correlations, 
including a diverging correlation length at the MSA critical point 0. 

However, it fails badly at low densities [|[ [|, and, as explained elsewhere H [|, [fl], the 
MSA (and GMSA) appears to be inferior to DH-based theories for describing the critical region. 
Hence, we have sought to remedy the deficiencies of the DH correlation functions as well, and 
specifically to do so in a more natural way. By following the spirit of the DH approximation, 
we extended the theory to the case of non-uniform, slowly varying ionic densities, p±(r), 
thus enabling derivation of a Helmholtz free energy functional M, [HJ. Ion correlations may 
then be obtained by functional differentiation techniques. This generalized Debye-Huckel 
(GDH) theory was applied to the calculation of Gnn(?), and provided not only the expected 
critical divergence of the second-moment density correlation length, £, but also the surprising, 
universal and exact divergence of £ in the low-density limit J8|, [l2| (where the GMSA fails 
by predicting a finite, non- universal value Q). This testament to the physical validity of the 
GDH approach motivated the calculation of the charge-charge correlations via GDH theory 
that is reported here. 

We find an expression for Gzz(r) = ( [p+(r) — p~(r)] [p+(0) — p-(0)] ) which in the low- 
density limit approaches the conventional and exact DH result, but which also explicitly satis- 
fies the Stillinger-Lovett second-moment condition. Furthermore, it exhibits charge oscillations 
for densities above a "Kirkwood line" in the density-temperature plane JtJ . More concretely, we 
find for the RPM (equisized hard-spheres with diameters a, charges iqo, and solvent dielectric 
constant D) the closed-form, Fourier transform expression for the charge-charge correlation 
function 

Gzz(k; Pn, T) = p N k 2 /[n 2 D + k 2 + a~ 2 gu(n D a, ka)], (1) 

where px = p+ + p- is the total ion density while the Debye parameter is given, as usual, by 
Kjj = 4irq 2 pN /DksT, and 

go(x, q) — x 2 (cosq - 1) - [21n(l + x) - 2x + x 2 ](cosq - sinq/q). (2) 

Expansion in powers of k yields Gzz(k; pn ,T) — (DkBT/Anq 2 ) k 2 + 0(k 4 ), which demon- 
strates satisfaction of both the Stillinger-Lovett second-moment condition, 



drr 2 G zz (r) = -6p N /Kh, (3) 
and the "zeroth-moment" or charge-neutrality condition, 



drG Z z(r) = => / drG ZZ {v) = -p Nl (4) 
J\r\>a 

for all pjv and T. In the low-density limit (|l|) becomes Gzz(k) ~ pNk 2 /(n 2 D + k 2 ), the 
exact, universal limiting behavior. By analyzing the poles of ([[]) we may obtain the predicted 
large-distance behavior of Gzz( r > PN>T): from that we find that simple exponential decay 
persists only up to x = n^a = xk] for x > xk the decay is oscillatory. Numerically we obtain 
the "Kirkwood value" 

x K ^ 1.17832, (5) 

which lies in the usually expected range ji], ||, ||, f?J . 

Before presenting the GDH calculation, we summarize briefly the conventional DH correla- 
tion functions for comparison. Debye-Huckel theory provides an approximate result for CT (r), 
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the mean electrostatic potential at r due to both a charge of type a fixed at the origin and 
the corresponding induced charge distribution 0] , namely, 

(j>v H (r) = q a /Dr - q a K D /D(l + K D a), r < a, 

= q a e K ° {r - a) /D{1 + k d o), r > a. (6) 

The DH thermodynamics uses only the self-potential lim r ^ H ( r ) ~~ Qa/Dr] — q a n^jD(\ + 
kdo). In a first effort to obtain correlations, one may supplement the treatment of Debye and 
Huckel with the formally exact statistical Poisson's equation, which for the RPM states [2(b)] 

V 2 ff (r) = -(4nq a /D)G zz (r)/p N , (a = ±). (7) 

This yields a G^p s with a simple screening decay that (i) satisfies (Q) for all pn and T, but 
(ii) violates the second-moment condition everywhere except in the low density, Koa —> 
limit. Furthermore, this Poisson route says nothing whatsoever about the density-density 
correlations. 

A second approach is to parallel the derivation of by putting (j>® H (r) into a Boltzmann 
form, g aT (r ± 0) = {p a (r)p T (0)) / p a p T ~ exp[-/3g CT T (r)] for r > a with g aT (r) = for r < a. 
Then both G™ z z and G™^ may be obtained, although the latter displays no sign of the proper 
critical behavior. However, a more glaring defect of this approach is that G^ z z (r) violates 
not only the second-moment condition, but also the charge-neutrality sum rule! [This follows 
readily from the inequality sinhec > x (when x > 0).] 

A third path, perhaps the most travelled in the literature, is to linearize the Boltzmann 
form to get g a r(^) — 1 — (3q a (j> T {r), for r > a. This gives the same charge correlations as does 
the Poisson route, but for G^ N the Coulombic terms cancel completely, as in the MSA. 

Now, although (Q) is an exact relation, approximate theories are generally inconsistent with 
respect to some identities; indeed, only the exact solution can satisfy all possible relations. 
Our GDH theory || provides an alternate, formally exact approach to the correlations, which 
more closely follows the thermodynamic theory and may then be judged on its relative merits. 

Our guiding motivation is that a free-energy-functional formulation ensures density cor- 
relations that are sensitive to the critical point (i.e., the compressibility relation is satisfied 
by construction); thus, we are led to develop a DH theory for an inhomogeneous electrolyte. 
As detailed elsewhere fj(| , the approach is sufficiently general to allow the calculation, via 
functional differentiation techniques, of all the a, r correlation functions G CTr (r) for an arbi- 
trary multi-component, equisized hard-sphere electrolyte. For brevity, however, we restrict 
consideration here to the RPM, for which the imposed density variations 

P±(r) =p±[l± Acosk-r], (8) 

(i.e., in response to an external potential) result in the Helmholtz free energy 

/Mr)] - f[p a ] = -ip^G^(k)A 2 + 0(A 4 ), (9) 

where / = —(3F/V ||, Thus Gzz may be found by expansion in A. 
The generalized Debye charging process yields the free energy [nol n3| 



F[ Pa {v)\ =F HC [ Pa {r)]+ dX^ J dr' Pa (r')^(r',{\q„}). 



(10) 



where F HC denotes the pure hard-core Helmholtz free energy functional while the mean 
electrostatic potential seen by an ion of type a fixed at r' is found from 

V CT (r') = Hm [<k(r;r') - q a /D\v - r'|], (11) 
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in which 4>„(r; i*') is the potential at r due to both the fixed charge at r' and the induced charge 
distribution [|[ |l(J. [Compare with (^) above, et seq.] To calculate 4> a (v; r') we begin with the 
exact inhomogeneous version of the statistical Poisson's equation (H}, 

V^ CT (r; r') = -(4tt/£») £ ? T/ o r (r)g TCT (r; r'). (12) 

T 

The DH approximation is to replace the g T(J with Boltzmann factors which depend on the 
potential. However, it is crucial to note that the varying ionic charge densities carry along an 
imposed overall electrostatic potential $(r), determined simply by 

V 2 $(r) = -4wq oPz (r)/D = -(Att/D) £ q T p T (r), (13) 

with appropriate boundary conditions. This is independent of the fixed charge of type a at 
r', and therefore should not contribute to the Boltzmann factor for g Ta . Hence, in the spirit 
of DH, we take 

g T<T (r; r') ~ exp[-/3g T ^(r; r')], |r - r'| > a, (14) 

with the "local induced potential" 

^(r;r') = ^(r;r')-$(r), (15) 

and, as before, g Ta = for |r — r'| < a. The need for the separation of 4> a (r;r') into a 
background $(r) and an induced piece <fi a (r;r') is clear in the limit |r — r'| — > oo, in which 
ln[g Tcr (r; r')] must vanish while CT (r;r') — > $(r) 

Next one inserts ( [Ti] ) into (|lj) and makes the second approximation of the DH procedure, 
namely linearization. This results in the full GDH equation 

V 2 ^(r;r') - -(4n/D) [g CT <5(r - r') - q p z (r)] , |r-r'|<a, 

= ^(r)^(r;r'), |r - r'| > a, (16) 

where the spatially varying Debye parameter H |) is defined by K 2 D (r) = (in/D) J2 T QtPt( v )- 
The second term on the righthand side for r — r'| < a, i.e., inside the hard sphere, is needed 
to cancel the background charge density pz( r ) there; it represents an effective "cavity source" 
term. 

To obtain Gzz for the RPM we chose the spatially varying densities (||), for which h 2 D {r) 
reduces simply to K 2 n . The resulting GDH equation (^) can be solved readily by Green's 
function methods ]^,]l0|. Integrating the self-potential Vv(r') against the density and charging 
according to (^) gives the free energy to order A 2 , from which our main result (^) follows by 
use of (||). Note that the q a S(r — r') source does not contribute directly to the charge-charge 
correlation function in the simple case of the RPM; rather it is the cavity term and <&(r) that 
serve to determine Gzz- 

To elucidate the long distance behavior of Gzz{y) we solve for the pole, ko, of Gzz (k) 
that lies closest to the origin in the complex k plane. The real and imaginary parts of kg 
plotted in fig. 1 were found by solving the coupled equations Re[G^(fco)] = ^ m [Gzz(^)] = 
numerically, using the Newton-Raphson method. When this pole is purely imaginary, 
corresponding to the leftmost part of curve (a) in fig. 1, Gzz decays monotonically as e~ r ^ z /r, 
where the screening length is £z — l/Im(fco)- In the low-density limit one finds £^ = 
«d[1 + \x 2 — ix 3 + + . . .] so that £z approaches the Debye value £d = 1/kd'- see 
curve (c). As kd and pn/T increase, the pole fco and a nearby, purely imaginary pole 
fci, curve (d), approach: at the Kirkwood value, n^a = xk [see (ph], they merge, with 
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Fig. 1. - Plots versus x = (Anq^pNa 2 / DksT) 1 ^ 2 of (a) Im(fcoa) = a/£z, the inverse charge- 
charge correlation length, (b) Re(fcoa) = 27ra/A, the inverse charge oscillation wavelength, 
(c) x = kdci, the inverse Debye length, and (d) Im(fcia), the subleading pole. Note that 
within pure DH theory the critical point occurs at x c = 1 and T* — 1/16; in pairing theories 
based on DH theory x c = 0.9-1.1 and T* = 0.052-0.057 [1,3,15]. 



a/^z — 2.266o, and for larger kd they become complex, implying the oscillatory decay 
Gzz{?) ~ cos[(27rr/A) — 9 KD ]e~ r ^ z /r, with A = 27r/Re(fc ) [see plot (b)] and 6 KD a phase 
shift. Hence, charge oscillations occur for densities p* N = pno? > x 2 K T* /Att ~ 0.110 T*, where 
T* = DkBTa/q 2 . Finally, at x = xx — 6.65232 the poles move to the real axis, i.e., they 
merge with their complex conjugates. Here the oscillations, with wavelength Ax ~ 1.847a, 
are no longer damped; this is suggestive of the onset of crystallization |(| although this region 
certainly lies beyond the limits of validity of our approximation. 

As a consequence of our analysis one sees that the GDH theory for ion correlations is 
the only one available that satisfies the requisite sum rules (||) and (Q), gives exact results 
for Gzz( r ) an d GjvAr(r) in the low density limit || [l^], and behaves sensibly in the critical 
region, predicting (mean-field-like) diverging density fluctuations j^]. 

Finally, we remark on the addition of dipolar ion pairing. Following [[| and one may 
straightforwardly add pairing to the calculation of Gzz- However, in this case, as opposed 
to the density-density correlations, one expects there to be no major contribution from the 
pairs, since, in the center-of-mass approximation, the dipolar ion pairs, of density p 2 , appear 
as neutral objects that cannot contribute directly to the charge-charge correlations. One 
important role of pairing, however, concerns the location of the Kirkwood line in the (pn,T) 
plane: owing to the depletion of the free ion density, pi = p + + p_ = p^ — 2p2, the onset 
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of charge oscillations is expected to lie close to the critical region: see |llj, |15| for the critical 
parameters and degree of pairing in the various approximations. Lastly, we mention that the 
charge-charge correlations in the GDH formulation with pairing still satisfy the sum rules (||) 
and (Q); however, in the Stillinger-Lovett rule (||) the "background" dielectric constant D 
that enters the definition of Kr> acquires a linear dependence on P2{pn ,T). In general, the 
state-dependence of D seems an open question |Q although most authors seem to hold that 
it should be totally absent in the RPM. If that is correct the pairing treatment would require 
further improvement. (See also |Tl}|.) 
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